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Abstract 



Under special conditions bacteria excrete an attractant and aggregate. 
The high density regions initiahy cohapse into cyhndrical structures, 
which subsequently destabilize and break up into spherical aggregates. 
This paper presents a theoretical description of the process, from the 
structure of the collapsing cylinder to the spacing of the final aggregates. 
We show that cylindrical collapse involves a delicate balance in which 
bacterial attraction and diffusion nearly cancel, leading to corrections to 
the collapse laws expected from dimensional analysis. The instability of 
a collapsing cylinder is composed of two distinct stages: Initially, slow 
modulations to the cylinder develop, which correspond to a variation 
of the collapse time along the cylinder axis. Ultimately, one point on 
the cylinder pinches off. At this final stage of the instability, a front 
propagates from the pinch into the remainder of the cylinder. The 
spacing of the resulting spherical aggregates is determined by the front 
propagation. 
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The formation of a singularity — the divergence of a physical quantity in finite 
time — is central to diverse fields including nonlinear optics, gravitational collapse, 
and fluid mechanics. The structure of singularities has been worked out in many 
examples for which a physical quantity blows up at a spatial point Typically, 
singular dynamics are self-similar: the characteristic scale separation between the 
singular and regular parts of the solution leads to the slaving of the spatial structure 
to the time dependence via scaling laws. The situation can be more complicated when 
many singularities form collectively and simultaneously. In this paper we analyze a 
simple example for which multiple singularities form in a short time. This work was 
motivated by a recent experiment in bacterial chemotaxis 1^-0. 

The experimental observation is shown in Figure |l]. In the first panel a diffuse 
cloud of Escherichia coli covers the depth of a petri dish filled with agar. The envi- 
ronment is prepared so that the E. coli excrete an attractant; each bacterium attracts 
all the other bacteria, and a cloud can collapse. In the second panel, the diffuse cloud 
collapses as a cylindrical structure, with highest bacterial density on the cylinder 
axis. In the final panel, the cylinder breaks down into spherical aggregates. In this 
paper we analyze this process, by constructing a similarity solution to describe the 
cylindrical collapse of bacteria, and then analyzing its stability. 
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FIGURES 




FIG. 1. Experiment showing formation and instability of a collapsing bacterial cylinder 
(Reproduced from Q ) . The first panel shows a diffuse cloud of bacteria filling the depth of 
a petri dish filled with agar, which then collapses (second panel) into a cylindrical structure. 
The cylinder subsequently destablizes into spherical aggregates. Details of the experiments 
are described in Budrene and Berg |^J6|. 



Chemotaxis in E. coli is an excellent model system for studying singularity for- 
mation. The biochemical response of E. coli to a changing environment has been 
thoroughly characterized IPHT^ over the past twenty-five years, so we understand 
how the bacteria sense and respond to their environment. As a consequence, it is 
possible to write down a "first principles" hydrodynamic theory for the motion of 
many bacteria [0,113 in which the response coefficients are measurable. Quantitative 



comparison between theory and experiment is possible, and any discrepancies can 
be traced directly to the biochemistry of individual bacteria [0. The application to 
singularity formation arose from the recent discovery by Budrene and Berg [^,^ of 
an assay in which E. coli excrete aspartate, an amino acid which is also an attractant 
for nearby bacteria. Attractant diffusion drives aggregation because it leads to an 
effective force between individual bacteria: a higher density of bacteria in a given 
region leads to a higher attractant concentration, which then attracts more bacteria. 

The initial interest in the Budrene-Berg experiment was stimulated by the sym- 
metrical patterns that form when chemotactic bacteria are seeded in the center of a 
petri dish, as shown in their papers 

Several theories have been developed for these patterns, most of which |T^-^ 
view the pattern formation as resulting from a linear instability of a (1 dimensional) 
travelling wave of bacteria. Recently, it was pointed out [0] that each of the aggregates 
in a a pattern corresponds to a density singularity in the hydrodynamic description 
of the bacteria. Therefore the pattern formation depends crucially on the dynamics 
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of singularity formation. Singularities in chemotaxis were anticipated by Nanudjiah 



^ and Childress and Percus p2| in studies of mathematical models of chemotaxis. 
An important feature, understood first by Childress and Percus, is that chemotactic 
collapse has a critical dimension: although collapse to an infinite density sheet is 
mathematically impossible, collapse to infinite density lines and points both can occur. 
It was argued in |^ that these facts crucially affect the patterns that can form . 

In particular. Figure 1 shows a step in the formation of the aggregates. The 
initially diffuse band (filling the depth of agar) cannot form a singularity by collaps- 
ing only one of its dimensions to zero thickness; instead it collapses into a cylinder 
(contracting two of its dimensions simultaneously). The cylinder later destabilizes 
to form aggregates, for which all three dimensions contract simultaneously. Models 
T5|^0[ viewing aggregate formation as the linear instability of a band cannot account 



for these experimental observations. These two different pictures of aggregates form 
lead to different conclusions about which biochemical parameters set the wavelength 
and structure of the patterns. For the "collapsing cylinder" mechanism advocated 
here, the characteristics of the pattern are set by the same biochemical cutoff which 
prevents an aggregate from reaching infinite density. 

Cylindrical collapse is also important when a uniform density cloud of bacteria 
breaks into aggregates. Linear stability analysis of the uniform density state predicts 
that the cloud directly breaks down into spherical aggregates. However, experiments 
p3| find that the the clumping is hierarchical: the uniform density cloud first collapses 
as cylindrical structures which then break into spherical aggregates. An important 
unsolved question is to explain the geometry of the high density regions during col- 
lapse, and to predict the distribution of final aggregates. 

In this paper, we use a combination of simulations and asymptotics to describe 
the breakdown of a cylinder in of three principal steps (Figure First, the bacteria 
collapse as a cylinder towards a line of infinite density. In the second step, uniformity 
along the cylinder axis is broken, and a singularity develops at a single point. Finally, 
the remaining cylinder breaks up producing a sequence of spherical aggregates. Our 
primary conclusion connecting the present theory with the experiments is that the 
spacing between the aggregates in patterns such as Figure 2 is determined by the 
local depletion of chemicals which make aspartate production possible. According to 
Budrene , the most likely candidate for this is the overhead oxygen concentration 



in the cell. This prediction is qualitatively in accord with the experiments; moreover 
this parameter dependence could be directly tested in future experiments, and would 
serve to discriminate this theory from those based on pure linear stability analysis. 
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FIG. 2. Schematic of the three phases of cyUndrical collapse. Top, a collapsing cylinder 
is formed. Middle, the cylinder becomes modulated. Bottom, the modulated cylinder 
pinches off and contracts, leaving a series of spherical aggregates. Note this sketch does not 
show the radial contraction of the cylinder which takes place as it collapses. 



In the following section, we review the basics of chemotaxis, and discuss details 
necessary to understand the Budrene-Berg experiments. Section 2 describes the cylin- 
drical collapse of the bacteria. Cylindrical collapse (in which two dimensions of the 
cloud contract simultaneously) is a critical case p2| , p^ , for which diffusion and at- 
traction nearly exactly balance. This criticality has complicated previous attempts 
to solve the collapse; a recent treatment of collapse in the 2D nonlinear Schrodinger 
equation inspired our solution. In Section 3 we perform a stability analysis of 
a collapsing cylinder. Perturbations to the cylinder can be described by a "phase 



equation" for the singular time [^. Solutions to this envelope equation explain full 
numerical simulations of a modulated cylinder. Section 4 describes the final stage 
of the breakdown of the cylinder, after the cylinder has pinched off at a point. Sta- 
bility analysis predicts the breakup of the remaining column of bacteria, a situation 
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analogous to a propagating Rayleigh instability in a liquid column ||2^. Appendices 
describe our numerical methods, and fill in the details of the calculations. 



I. PROBLEM FORMULATION: BACTERIAL CHEMOTAXIS 



Chemotaxis refers to the migration of bacteria up chemical gradients. For Es- 



cherichia coli, the basis of chemotaxis is largely understood P8| , p9| : in the absence 
of a chemical gradient, an E. coli bacterium performs a random walk When 
chemical gradients are present, the bacterium's internal biochemical reactions detect 
the gradients and couple to the bacterial movement system. This sensing biases the 
random walk, and the bacterium has a net drift towards a chemical attractant. Un- 
der special conditions, the bacteria can excrete the chemoattractant aspartate 
by converting carbon and nitrogen sources (succinate and ammonia, respectively) in 
their environment. In these experiments no external chemical gradients are present. 
Instead, each bacterium moves in response to the attractant produced by other bac- 
teria. Thus, the excretion of attractant produces a long-range force between the 
bacteria, and induces complicated interactions in the colony. 

The equations for the collective motion of the bacteria can be derived (with no free 
parameters) from the underlying biochemistry allowing quantitative comparisons 



between theory and experiments. The basic equations for the bacterial density p and 
the attractant concentration c are 

9P r^^2. 



= DbV'p - V • [kpVcj + ap (1) 
^ = D^W'c + ap. (2) 

Here Di, is the bacterial diffusion constant, k the chemotactic coefficient, a the rate of 
bacterial division, a the rate of attractant production, and the chemical diffusion 
constant. The terms in equation ^ include the diffusion of bacteria, chemotactactic 
drift and division of bacteria. Equation |] expresses the diffusion and production of 
attractant. 

Equations of this type were first used to describe bacteria by Keller and Segal 
and, with variations, have been the subject of extensive investigations (see, e.g. 



3^ , |33[] ). For E. coli, Schnitzer et. al. established the connection [0,113 between the 
time averaged properties of the bacterial response and parameters in equations (|l], |^). 
Thus, the extensive studies of individual bacteria provide a rigorous justification for 
the equations, as well as measurements of the coefficients. There is one complication 
to this statement that is worthwhile to mention: in 1975 Spudich and Koshland 



showed that E. coli have "non-genetic individuality", manifest in a distribution of 
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tumble times (by about a factor of 2) between genetically identical bacteria. The 
consequence of this is that both time-averaged and ensemble-averaged properties of 
the bacteria are necessary to predict hydrodynamic coefficients; in addition, the dy- 
namics must be such that the "distribution" of bacteria in the ensemble does not 
change with time. 

It is convenient to nondimensionalize equations (||, ^ by choosing a character- 
istic density equal to the maximum initial density po- The characteristic scale of 
attractant is Dt/k. The density then determines the length scale and timescale ac- 
cording to H = DhDc/ {akpo) and to = Dd (akpo). Typical numerical values are Dj, = 
7xl0~^cm^/sec; D^. = 10~^cm^/sec; k = 10~^^cm^/sec; and a = 10^/second/bacteria. 
For an experiment which has po = 10^/ cm^, the length scale is 260 microns and 
the timescale 100 seconds. The equations become 

1^ = VV - V ■ (pVc) + 6p (3) 
^l = V^c + p, (4) 



where e = Db/Dc and 6 = ato- For the experiments shown in Figure 1, cells divide 
much more slowly than the dynamics take place: the timescale for cell division is ~ 2 
hours, giving 6 ~ 0.01. Therefore we approximate 6 = 0. The value of the parameter 
e varies: for experiments in semi-solid agar, the diffusion of bacteria is much slower 
than attractant diffusion, which motivates the limit e = |^. For experiments on 
bacteria in a liquid culture, e ~ 1. We will consider both of these limits in this 
paper. The e = limit is particularly convenient for asymptotic calculations; we will 
use e = when doing analytic calculations. Our numerical simulations give results 
independent of e in the range between and 1. 

For analytic calculations, working with the mass is particularly convenient, as 
will become apparent below. For reference, we show the form of the equations here. 
Define the mass contained within a radius r as 

m(r) = dr r'^~^p. 



This definition (and the limit e = 0) allows us to eliminate the concentration and 
write the original equations as 



7 



II. CYLINDRICAL COLLAPSE 



In this section we construct the solution for a uniformly collapsing cylinder of 
bacteria, according to the evolution equations (3,4). Although the existence of cylin- 
drical chemotactic collapse is well known p4| , the asymptotic solution describing the 
collapse has never been constructed. Herrero and Velazquez p3 understood impor- 



tant features of two-dimensional collapse and attempted to find the solution; however, 
their solution is different from what we observe in numerical simulations (see Figure 
^ and is therefore probably unstable. 

The usual derivation of a similarity solution derives a set of ordinary differen- 
tial equations from the scaling laws suggested by dimensional analysis. As shown 
below, the major difficulty arises because in the present situation these ODEs have 
no solutions consistent with the boundary conditions. This breakdown of "dimen- 
sional" scaling only occurs in two-dimensional collapse. In higher dimensions, the 
scaling laws suggested by dimensional analysis work fine The overall structure 
of this problem is similar to critical collapse in the focusing nonlinear Schrodinger 
equation where the initial proposal for constructing the asymptotic solution in- 
volved rather intricate mathematics (continuation as a function of spatial dimension). 



Recently, a more physical derivation was put forward |2^. This work on the non- 
linear Schrodinger equation inspired our construction of the solution for chemotactic 
collapse. Our central result is that the density singularity shows corrections to the 
dimensional scaling laws. Defining t = t* — t the distance to the singular time, the 
dimensional scaling law is p = r~^. In the transient regime, we find a slow power law 
correction to this dimensional scaling law, p ~ r~^^/^. After the transient asymptotic 
regime, logarithmic corrections are present, of the form 



P 



log|logr|| 



r 



Although the transient regime is present for 10 decades in the collapse time, we believe 
that this log log law is the asymptotically correct blow up rate. A log log correction 



also appears for 2D critical collapse in the nonlinear Schrodinger equation (See p5 
and the references therein). 



A. Critical Dimension for Collapse 



The competition between dissipation and collapse leads to a critical dimension. 
In this problem, the critical dimension is 2, and one dimensional collapse — that is, 
collapse to a planar structure with infinite density — is forbidden. 
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We make qualitative arguments to explain the critical dimension by comparing 
the chemotactic and diffusive ffuxes in a contracting structure. First, no singularities 
exist for one-dimensional contraction. For a sheet of thickness C. the diffusive flux is 
of order (see equations ([l|, |^)) 

~ -d/j. (6) 

The chemotactic flux follows by integrating Dcc" ^ ap and defining M^^ as the mass 
per unit area of the planar region. Then the chemotactic flux is 

^ fcpc' ~ akpM^^D;\ (7) 



If the system collapses onto a plane, the thickness of the sheet £ — > 0. In this case the 
diffusive flux blows up while the chemotactic flux is unchanged. Thus a planar region 
with small thickness is unable to form infinite density, because diffusion eventually 
stops the collapse. 

The situation is different for higher dimensional structures. For symmetric spher- 
ical collapse (three directions contract simultaneously), the chemotactic flux is sin- 
gular. When we balance DcV^c ~ ap, we find (r^c')' ~ ar'^p/Dc- This implies a 
concentration gradient c' ~ aM^^ /{i'^Dc), where M^^ is the mass contained within 
a radius i. The net inward flux of bacteria is then 

-Dbp akpM^^D-^ 



As £ — > the inward flux (second term) dominates and collapse occurs. 

In two dimensions we encounter a subtlety. Assuming cylindrical collapse and 
repeating the dimensional argument, we have (rc')' ~ arp/Dc, and c' ~ aM"^^ / {iDc). 
(Here M^^ is the mass per unit length of the cylinder.) The inward flux is 

_ -DbP + akpM^^p-^ 



Two dimensional collapse is critical: both fluxes have the same scaling with l. Ac- 
cording to this simplified argument, there is a net inward flux if M"^^ > DcDi)/{ak), 
which suggests that a system with mass above this critical value collapses. 



B. Similarity Solutions 



We now quantify the preceding dimensional arguments and collect the known so- 
lutions to the chemotactic equations. (As discussed above, the analytic solutions are 
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derived with e = 0.) First, consider one- dimensional collapse. Making the substitu- 
tion ^ V = Vc = dxC in equations H) implies 

dv d'^v dv 

— V- 



dt dx^ dx 

This equation is the Burgers' equation; singular solutions to this equation do not exist 

In d = 2 and higher, density singularities can develop. In three dimensions, the 
nature of the blowup is straightforward. The characteristic length scale (L) varies 
in time, and the spatial structure is determined by the changes in L. A singularity 
corresponds to L ^ 0. We guess the form of the similarity solution by balancing 
the different terms in equations The diffusive dynamics imply L = y/t* —t = 

^/T, with t* the singular time and r the distance to the singular time. Defining a 
dimensionless similarity variable rj = r/L, we find the scaling form of the density, 
concentration, and mass: 

p = ^Riv) 

c = Civ) 
m = L'^-^M{7]) 

In writing this form of solution, we have assumed radial collapse at the origin (r = 0). 
For a similarity solution to be valid, it must obey the correct boundary conditions: 
the density p and the attractant concentration c must be time-independent far from 
the singularity, which requires R ~ r^"^ and C ~ constant as rj ^ oo. 

Plugging in the scaling form gives an ordinary differential equation in the similarity 
variable rj. (Throughout this discussion, d refers to the number of simultaneously 
contracting dimensions.) 

uML = ^d-^R' + RM (8) 

r]'^-^R = M'. (9) 

In (i = 3, the similarity equation can be solved exactly; the one stable solution, found 
by Kadanoff [p6| , 



IS 



4(3 + r/2^ 



R-^r^. (10) 



As demonstrated in [RUl, this solution well describes numerical solutions. 
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FIG. 3. Density profiles for a collapsing cylinder. Different lines correspond to different 
times approaching the singularity. On the right, the density has been rescaled by the 
maximum density and the radius rescaled by the width. The solid line shows a profile at 
the same time both in this Figure and in Figure ^. The crossover between the stationary 
solution, for which the large rj behavior scales as r/~^ and the true outer solution, which 
scales is r]~^, is apparent. 



C. Solution in two dimensions 



In two dimensions, there is no solution to equations (|,^ which satisfies the bound- 
ary conditions. To see this, note that the similarity equations can be integrated to 
give 

R = e-e-i~ (11) 

This form for R cannot satisfy the boundary conditions that the density and mass be 
stationary at large rj, because R grows without bound a.s r] —>■ oo. 

Nevertheless, a similarity-type solution to the equations exists, as Figure || shows. 
The density profiles at different times (left) have been collapsed by rescaling the length 
scale (right). This figure shows the basic features of the solution: in the "inner" 
region, there is a scaling solution. This matches onto an "outer" region far from the 
singularity, with a final "boundary" region where the boundary conditions must be 
satisfied. We illustrate this section with plots from a single numerical simulation. 
The numerical method is described in Appendix B; its most important feature is the 
remeshing, which moves mesh points every time the maximum density increases by 
one percent to better resolve the singularity |^ . 
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How do we reconcile the numerical (and experimental) observations of two- 
dimensional collapse with the argument above that no solutions exist? We argue 
here that corrections to the similarity solution arise to solve this problem. We find 
in numerical simulations that although the basic self-similar scaling pm ~ is 
preserved, the time dependence of L is different than the simple dimensional argu- 
ment suggests. Figure ^ shows the scaling of L with time: if the dimensional law 
L ~ ^/T held, then L/y/r (plotted on the right) would be constant. Our simulations 
indicate a more complex time dependence than L ~ s/t. We argue that the prob- 
lem develops two timescales: a fast timescale on which collapse happens, and a slow 
timescale corresponding to changes in the similarity profile. Working in "similarity 
variables" allows a separation of these two timescales: we eliminate the known, fast 
time dependence so we can analyze the slow evolution of the system. 
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FIG. 4. Time dependence of the length scale, L. In simulations, the length scale is 
defined as the length over which the density decreases by a factor of five, a characteristic of 
the inner region. The first plot shows L versus distance to the singular time, which obeys 
close to a square root law. The second plot shows versus r, showing evidence of the 

transient regime [Lj yfr ~ r^^"^^) and, closer to the singularity, the asymptotic regime with 
logarithmic corrections to the length scale. For comparison, we have plotted the predictions 



of Herrero and Velazquez |35|. 
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As L — 0, the inner collapsing region converges to a similarity solution which we 
call the stationary solution because it has the same form as the stationary solution 
to the original equations. This asymptotic solution has a dimensionless mass of 
precisely 4. Thus, the evolution of the similarity solution has a specific physical 
interpretation. In dimensionless (similarity) variables, an inner region expels any 
excess mass to approach M = 4. (Figure Matching between the inner and outer 
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regions determines the dynamics. 
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FIG. 5. Schematic of cyhndrical cohapse. 



1. Separation of Time Scales and Inner Solution 



To find the solution, we use techniques originally applied to a similar problem in 
the 2D nonlinear Schrodinger equation |2^. We first define an important quantity: 



the collapse rate A{t) = —LL. (Note that the collapse rate of the system is L/L; 
we can make a dimensionless collapse rate by multiplying by the timescale t L'^. 
Thus A is the collapse rate.) For "regular" collapse, the collapse rate —LL = 1/2 is 
a constant. In the presence of corrections, the collapse rate goes asymptotically to 
zero. 

We rewrite the similarity equation assuming that both quantities M and R depend 
on the similarity variable rj and (slowly) on time. The similarity equation is then 

L'^^ + AriM' = rjR' + RM (12) 
dt 

riR = M'. (13) 



Note that here the time derivative of M refers only to explicit time dependence of the 
mass; the second term in the equation takes into account the time dependence slaved 
to the varying length scale. 

We can solve the similarity equation in the inner region. As the collapse proceeds, 
it does indeed slow down. In fact, a numerical measurement of A (Figure |^) shows 
that the collapse rate decreases by 3 orders of magnitude during the initial stages of 
the collapse. This motivates us to look for a solution with A = 0; that is, a stationary 
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solution (in similarity variables). This stationary solution solves Equations (|T2],|T3|) 
when A = and dtM = 0. The equation for the mass is then 

r]M" + M'{M - 1) = 0. 
Exact analytic solutions to this equation are 
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FIG. 6. Plot of the collapse rate A and r/^^ versus distance to the singular time r. The 
hne shows a power law A ~ r^^^. A is measured by evaluating LL numerically; L is the 
difference in length scale between two time steps, divided by the time step size. The length 
scale L is defined as the radial coordinate where the density decreases by a factor of 5. 
Similarly, r/* is measured as the coordinate where the numerical density profile differs from 
the exact stationary solution (Equation (p^)) by a factor of 4. 

From the formula for Mg, we see an important feature of the stationary solution: 
the total mass is 4 in dimensionless units. This reflects a rigorous result p4| , ^ : 
collapse will occur if and only if the total mass per unit length of the cylinder sat- 
isfies M > 4. For M < 4, no collapse is possible — the system evolves to constant 
density. For M > 4, collapse occurs; however, numerical simulations show that the 
collapse converges toward a solution with a collapsing mass precisely equal to four. 
In similarity variables, therefore, mass flows away from the origin. 
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Compare this expression for Ro to the numerical density profiles in Figure ^ The 
shape of the profile is accurately described by the formula for i?o, which confirms that 
the stationary solution holds in the inner region. At large rj, the stationary solution 
has Ro ~ 8?7~^, while the boundary conditions require R ~ cri~'^. This "inner" 
solution to the equations must therefore match onto an "outer" solution, as shown in 
Figure ^. How to analyze this matching is the subject of the rest of this section. 

Where does the crossover between inner and outer region occur? We can estimate 
the location of the crossover by noting that the crossover between inner and outer 



solutions will happen at some coordinate 77=,,, when ArjM' ~ rjR' (see Equation (JT^) 
This gives 

V. ~ A~'/\ (16) 
The proportionality between A and r/"^ is clear in numerics (Figure 0). 



2. Transient Regime 



Near the beginning of the collapse, the singular (inner) part of the solution has 
formed, but still strongly feels the influence of the boundary. During this initial, 
transient regime in the dynamics, excess mass is expelled. At the end of the transient 
regime, the true asymptotic state is reached. The transient regime is clear in the 
simulation shown here (Figures ^, ^ |^), where it persists until r ^ 10~^°. Here we 
attempt to analyze the transient regime and consider how the system approaches the 
asymptotic regime. This analysis is difficult, because the transient dynamics depend 
on the initial and boundary conditions of the system, and hence is particularly difficult 
to characterize analytically. 

Why even bother with the transient behavior? Because, as we mentioned above, 
the transient regime persists until r ^ lO^^*'. The density increases by 10 orders 



of magnitude before the true asymptotic dynamics are reached ||3^. Therefore, all 
experiments can probe only the transient regime. Furthermore, the detailed dynamics 
of the transient regime are not highly sensitive to the initial conditions. We did a 
number of simulations in which the length scale of the initial density distribution was 
varied (Figure |^). When the density is initially almost constant, the time at which 
the singularity is reached, t*, is large compared to a simulation in which the density 
is initially peaked. However, the transient dynamics are almost identical, regardless 
of whether the density is initially constant, or highly peaked [OT. 
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FIG. 7. Comparison of different initial conditions. On the left, the time until the singu- 
larity is reached depends logarithmically on the length scale of initial density variation, as 
expected for an unstable collapsing mode. On the right, the transient regime dynamics are 
almost identical, despite the wide range in length scale of perturbation. The legend shows 
the length scale of initial density variation for each curve. 



Our discussion of the transient regime that follows is incomplete and ad hoc. It 
is, however, the only self consistent explanation of numerical results we were able to 
find. We welcome additional, more conclusive work on this question. 

We argue that the corrections to the length scale in the transient regime are 
faster than logarithmic. Although theoretically we cannot rule out logarithms, we 
were not able to find an asymptotic solution for a correction which depended purely 
logaritmically in time that explained all features of the numerical simulations. For 
this reason we believe that the most consistent interpretation of the transient regime 
is that it actually breaks the scaling laws predicted by dimensional analysis and 
introduces different exponents. 

Here we use features of the numerical solution to guide our construction of the 
transient correction; this argument, while not entirely satisfactory, is at least consis- 
tent with all the simulations. We seek corrections to dimensional scaling of power 
law form: that is, L ~ \fTT^ and we are determining (5. As long as /3 > , we still 
have asymptotic validity of the analysis, because A goes to 0. (Note that /3 depends 
on the initial and boundary conditions, and is not universal.) 

We can distinguish between the transient and asymptotic regimes by looking at 
Figure This plot shows profiles of the normalized density (-R(O) = 1) times r/^. At 
later times, the curve shows a fiat region, demonstrating that R ~ r^"^ over a large 
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region. This is the hallmark of the true asymptotic regime: the stationary solution 
is valid in the collapsing region, then r^"^ outer region exists far from the singularity, 
which in turn connects to the boundary. For earlier times, however, the flattened 
region is not present: instead the forming singularity connects directly to the outer 
"bump" where the excess mass has been pushed. ETI 




Similarity Variable r\ 

FIG. 8. A plot of Rrj^ versus ry, used to construct the corrections in the transient regime. 
Profiles are for different times. The line shows a power law r]^R ~ r/^/^. The heavy line 
shows the profile at the same time as the heavy line in Figure ^. 



The inner solution must match on to the "bump" that makes up the outer region 
4^ . The bump corresponds to the excess mass in the system, which is stationary in 
real coordinates. Thus the position of the bump in similarity variables is rjb ~ ~ 
^-1/2-/3^ where we have assumed power law corrections to the length scale. By reading 
off the plot (Figure ^) we can see two important features of the outer region: first, 
the value of rj'^R is constant at the position of the bump: rj^Rt = 7- Next, the form 
of the solution from the crossover point to the bump [r]^ < 7] < rjh) is approximately 
a power law, which we write 

r]^R = of for 7]^, < r] < r]b 



The exponent a is a free parameter that we find from the simulations. Taking a from 
the numerics (reading off Figure |^) uniquely determines all other scaling exponents 
and hence gives a consistency check. To determine (3, combining the relations above 
lets us determine the time dependence c ~ ~ -j-o'i'^/'^+i^) _ This means that in the 
intermediate regime r]^ < r] < rjb the solution obeys 
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Finally, we match the two solutions at the crossover rj^,, leading to the result r/^, ~ 
^-a(i/2+/3)/(2+a)_ demanding consistency in time scaling we must have r]~'^ ~ yl ~ r^^. 
This determines f3: 



For the simulations shown, a = 3/4. We have examined a for a variety of initial 
conditions and find it always to be close to this value. This one assumption for 
a leads to predictions for other measureable exponents, all of which are in good 
agreement with numerical simulations. 



Thus in the transient regime we have 

A = r3/8 




L X 

FIG. 9. The scaling of the maximum density with time. On the left, pm versus L, which 
is well described by the expected scaling. On the right, pmT versus r, showing the transient 
and asymptotic regimes. 



3. Asymptotic Regime 



Once a true scale separation exists between the inner (collapsing) region and 
the boundary, the matching is different. Here we see slow (logarithmic) corrections in 
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numerics, motivating us to seek a logarithmic correction to the length scale. We define 
the slow time scale s = — logr, and seek a length scale L of the form L = ^/t/ f{s). 
With this modification, Equations ( |T^JT3| ) for the mass becomes 

1 dM 



f{s? ds 



+ ArjM' = r]R' + RM 



(17) 



The inner region must match onto an "outer" solution of the form R = cri~^ The value 
of the coefficient c can be estimated from the total mass and size of the system. If c 
is independent of time, the total mass M in the far field is logarithmically divergent: 
M ~ clogrj. However, conservation of mass requires that the total mass be constant 
at the boundary: M{W) = Mt, which implies c ~ Mt/ \og{W/L) (up to logarithmic 
corrections) or 

Mt 

c ~ . 

hgW + s 



As explained above, the crossover between inner and outer solutions occurs at rj^ ~ 
y4~^/^. At ?7* we require that the mass fiux dM/ds is continuous; this condition 
connects the inner and other solutions, and fixes the dynamics of /. 

From the similarity equation (ITTp and R{f]) = cri~'^ we have the mass fiux in the 
outer region 

^ = -Afc + 0{r^-') ry^oo (18) 



The mass fiux as r/ ^ oo is spatially uniform. This must match the mass expelled 
by the inner region, which we find from Mo ~ 4 — ?7~^. Equating fiuxes at rj^ gives 
dMo/ds\r„ = -Ape or 



We can convert this to an equation for / by using the definition oi A = —LL = 
/~^(l/2+/'/ /) ^ /~^/2, where a prime denotes derivative with respect to s. Plugging 
in the value of c from above, we have an equation of the form 

-f ^ Mt 
P \ogW + s 



This relation holds for large s, giving the scaling of / with time: 

-/' 1 

P s 
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which gives 

1 

''log I logr| 



/ 



The asymptotic solution for two-dimensional collapse is a very slow (log log r) correc- 
tion to the scaling laws from dimensional analysis. Our simulations cannot verify the 
functional form of /; however, we certainly see a crossover to slow (or nonexistent) 
corrections to the length scale. In Figures |^,^, and ^, the slow correction is visible 
for r between 10^^'^ and 10^^*^. A blow-up of this slow correction show that it is ap- 
proximately constant, but dominated by numerical noise — our numerics aren't good 
enough to resolve this correction. Note that a simple logarithmic correction (that is, 
/ proportional to log r) is not consistent with our simulation, because such a function 
would not be approximately constant over 10 orders of magnitude in time. Therefore, 
our numerics are consistent with, but not proof of, a log log correction. 



III. EVOLUTION OF A MODULATED CYLINDER 



A collapsing cylinder eventually breaks into spherical aggregates. Here, we find 
the envelope equation which describes how modulations to the cylinder evolve. The 
challenge is to describe a collapsing cylinder. Collapse amplifies initially small per- 
turbations; therefore small variations along the cylinder (in density and the radial 
length scale) become large. We can perform a valid perturbation analysis by study- 
ing variations in the singular time t*. 

In the original similarity solution, the singular time t* is undetermined: if t* 
changes by a constant, the solution remains valid. Allowing slow spatial variation 
in t* breaks this symmetry. Therefore we expect the variation of t* to produce slow 
dynamics in space and time. Because this mode is the most slowly decaying, it 
dominates the evolution of a cylinder. We derive a phase equation, an approach 
used in many problems when stability is governed by a slow mode associated with a 



broken symmetry [EQ]. Phase equations E3[ were invented to understand problems 



like convection, where the relevant symmetry is translation, and the stability analysis 
is relative to a travelling wave solution. Earlier research moved towards applying 
phase equations to singularities: In previous work on blowup in the semilinear heat 
equation, Keller and Lowengrub derived a transformation from a blowing-up 



variable to one that vanishes; they then perturbed in the vanishing variable. Also 
in the context of the semilinear heat equation, Bernoff MB has looked at how the 
singular time varies along a cylinder. 

We compare solutions of the phase equation to full numerics and show that the 
evolution of a modulated, collapsing cylinder is well described by the simplified phase 



20 



equation. 



A. Derivation of a Phase Equation 



This section gives a flavor of tlie derivation of tlie pliase equation, for details see 
Appendix 0. We seek an equation for the dynamics of t{z, t) = t* — t = to + T{z, t) —t. 
Slow variation requires that both T" <^ T' and T <^1. Once the collapse time varies 
along the axis of the cylinder, R and C are no longer exact solutions to the equations 
(0, 1). We write the correction as 

R = Ro{7]) + 5Ri{ri.z,t) 
C = Co{v) + SCi{v,z,t) 



where 77 = r/Lr, with the radial length scale Lj. = ^T{z,t)/ /{t), where / is the 
correction to the length scale. The perturbation parameter 6 is of order Lr/L^, where 
Lz is the scale of the density variation along the axis of the cylinder. 

We insert this guess into the original equations, and expand in 6. The lowest-order 
equation gives the original similarity equations. At first order, we find an equation of 
the form: 

A(i?i, Ci) = F{Ro, Co)iTt + 1) + GiRo, Co)tz. + H{Ro, Co)- (20) 



On the left side is a linear operator A acting on Ri and Ci; A comes from the 
linearization of the original equations. The right hand side contains derivatives of the 
singular time multiplied by known functions of Ro and Co- 

Although Rl and Ci are unknown, the right hand side is constrained by a solv- 
ability condition: Any function which is annihilated by the adjoint of A must be 
orthogonal to the right hand side. That is, if N^g = for a nonzero g, then the inner 
product 1^ 



{g,A{Ri,Ci)) = {A^g,{Ri,Ci)) = 0. 



Hence, the right hand side of (|20|) is orthogonal to g. In this case (see Appendix 0) 
precisely one nonzero g satisfies A"!" = 0. Taking the inner products leads to a phase 
equation of the form 

Ci(n+ l)+C2nz + C3— = 
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where the constants Ci, C2, C3 can be expressed in terms of the known functions g, F, 
G, and H: 



ci = {g,F{R,,Co)) 
C2 = {g,G{Ro,Co)) 
cs = {g,H{Ro,Co)) 



As shown in appendix Q for a cyhnder disintegrating into spheres the phase equation 
is 

—75— =^zz (21) 



The correction terms in this equation arise directly from the corresponding corrections 
to the length scale in the similarity solution. In the subsequent section, we show that 
this results in asymptotically different scalings for the radial and axial length scales 
on the collapsing cylinder. (Hence, a "point" singularitity that forms on a collapsing 
cylinder does not have the same collapse rate as spherical collapse.) In the absence 
of corrections to dimensional scaling, the phase equation is simply 

n + 1 = r,, - ^ (22) 



B. Numerical Simulations of a Collapsing Cylinder 



Now we compare solutions to equation (^) with a fully nonlinear simulation of a 
collapsing cylinder. We have found two different mechanisms by which modulations of 
the cylinder can produce singularities: The first is a "point" singularity, in which the 
density blows up at a point on the cylinder; the second is a "travelling" singularity, 
which moves along the cylinder axis with a diverging velocity (as the singularity is 
reached.) 

The primary technical difficulty in simulating a collapsing cylinder is developing 
a remeshing algorithm to closely approach the singularity. The meshing algorithm 
described here resolves density singularities along the axis of the cylinder with es- 
sentially arbitrary resolution. The algorithm is based on a simple one-dimensional 
scheme, which redistributes mesh points every 50 timesteps to resolve the singu- 
larity. The two-dimensional algorithm uses the one-dimensional remeshing scheme 
along both f and z simultaneously; the two dimensional equations are then solved by 
operator splitting. Details of the algorithm are summarized in Appendix 0. 



22 



A typical simulation |^ started with a z-independent initial condition, which 
was allowed to progress until the maximum density reached IC^. At this point, the 
radial profile of the collapse was well approximated by the (2-dimensional) collapsing 
similarity solution constructed in the previous section. We then added a 2;-dependent 
perturbation to the density profile, with amplitude much smaller than the ambient 
density. 

The separation of scales hypothesis underlying the derivation of the phase equation 
is maintained uniformly in time. We experimented with different functional forms for 
the density perturbations; as long as the length scale for variation in the z direction 
is longer than the variation in the radial direction L,,, perturbations tended to grow. 
In all cases, the relation ^ was maintained. As demonstrated in Figure |10|, a 



radial cross section of the cylinder always revealed density profiles in agreement with 
the two dimensional collapsing singular solution constructed in the previous section. 




FIG. 10. Figure showing the radial density profile at a single point on a modulated, 
collapsing cylinder. Note the region matched onto the region: The solution is 
well-described by the two-dimensional similarity solution throughout the collapse (cf. Figure 
0) even in the presence of a z dependent modulation. 



Extensive simulations have found two types of density singularities caused by 
modulations to the cylinder. 



1. Travelling Singularity 



Travelling singularities occur when a step-like perturbation is placed on the cylin- 
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der, increasing the density for z > zq and decreasing the density for z < zq. The 
subsequent evolution occurs at the boundary between these two regions. A simula- 
tion of this process is shown in Figure |1T[ The higher density region propagates to 
the left. Heuristically, the higher density region is beginning to contract as a sphere, 
so its decrease in size is consistent with the beginnings of spherical collapse. 

This propagating singularity can be described as a solution to the phase equation 
of the form 

T = To (i){z - Zo{t)) 

where Tq = t* — t is the basic phase expected from collapse. All the nontrivial space 
and time dependence is absorbed in cj) and Zo- Inserting into the phase equation (|2T|), 
we have 

— p — "Tj 

The right hand side of this equation is 0(ro), which becomes arbitrarily small as the 
singularity is approached. Therefore, the left hand side must be equal to zero, which 
gives 

fo0 - ZoToCj)' +1 = 0. 

If we demand the balance ZqTq = —A, then 

Zo = AlogTQ. 

The solution for (p is then 

0(r^) = 1 + e'?/^. 
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FIG. 11. Time evolution of the centerline density p{r = 0, z), with a step function initial 
condition. The perturbation is seeded when p = 10^, and subsequent profiles show when 
the maximum density increases by a factor of five. The highest density region propagates 
to the left. Hence, the spatial extent of the high density region shrinks as the collapse is 
approached. 



Figure |TT] shows the density along the centerline of the cylinder. The decay of the 
highest density is exponential, as predicted by the phase-equation solution constructed 
above. A fit to the numerical data shows that 



p(r = 0, z) 



Figure 12 shows the location of the edge of the maximum density region zo{t) as 
a function of r. As predicted by the phase equation analysis, the figure shows that 
Zq = AlogT. A least squares regression gives the prefactor A ^ 0.17. The qualitative 
features of the numerical simulations are thus in good agreement with the theory. 
Quantitatively, however there is a discrepency: the theory predicts that we should 
have A/zq = 1, while our numerical simulation gives A/zq ~ 0.6. We believe that the 
discrepancy arises because the convergence of our numerical scheme requires we keep 
a small e ~ 0.1, while the analysis in the previous subsection assumes e = 0. 
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FIG. 12. Location of the edge of the maximum density region zq(t) as a function of 
time. As predicted by the theory, the edge moves according to the law zq^t) = ^logr. 
Regression gives A ~ 0.17. 



2. Point Singularity 



Stationary singularities — in which the blow-up happens at a spatial point — occur 
in the numerics for a wide variety of initial conditions. We believe that this represents 
the generic evolution of a modulated cylinder. Indeed, it is the end state of the 
travelling solution just discussed, when the propagating wave runs into the reflection- 



symmetric boundary condition. An instance of this solution is depicted in figure 
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FIG. 13. Time evolution of the centerline density p{r = 0,z), for a point singularity. 
The singularity was initiated by placing a perturbation (symmetric about z = 10) on a 
uniformly collapsing cylindrical solution. 

In contrast to the previous case, the corrections to dimensional scaling are impor- 
tant in this solution. For our analysis, we take the point of blow-up to be at 2; = 0. 
Satisfying the equation requires that the z length scale incorporate corrections to the 
scaling: 



r = ToC, 



To^h{T) 



Using ^ = z/{t2) the phase equation becomes 

':-2 ( , 



r 



Demanding that the two sides scale the same way in time (and assuming h has a 
power law form, so that Toh' /h =constant), we have 



1 

P 



T, 



1-27 



/i2 



which gives 7 = 1/2 and f = h. Thus ~ t^^To which differs from the radial 
length scale L,. ~ t^^To''^^. The result is 



_-3/8 



(23) 
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This shows that the generic density singularity that forms during the breakdown of 
cyhndrical collapse, is not spherical collapse, but something milder. Locally, since 
the axial scale Lz is much larger than the radial scale Lr the structure still looks like 
a cylinder. Numerical evidence for this conclusion is evident in Fig 13, which shows 
that the singularity develops a length scale in the axial direction which is much larger 
than the radial scale 1/a/p. For example, in Fig 13 when p = 10^° (so Lr = 10~^) 
then Lz ~ 10~^. Our numerical algorithms have unfortunately not allowed us to find 
the asymptotic Lz/L^ numerically; the problem is that the separation of scales is so 
great between the radial and axial scales that one needs many more mesh points than 
we can afford to resolve the asymptotic regime. 



IV. BREAKUP INTO SPHERICAL AGGREGATES 



The question of relevance to the experiments is what happens next: Once blowup 
occurs at a spatial point the cylinder has a "free end", which changes the nature of 
the collapse. We can no longer use the strategy of the previous section — perturbation 
about a collapsing cyhnder — because the radial structure is no longer closely approxi- 
mated by the cylindrical solution. The pinch-off drives the dynamics; specifically, the 
sharp end of the cylinder forms a travelling wave. Heuristically, note that an edge of 
bacteria produces a higher concentration of attractant where the density is higher. 
Thus the "tail" of bacteria moves toward higher attractant density, and a travelling 
wave can form. Quantitatively, recall that for variation in one spatial dimension 
the equations reduce to the Burgers' equation, which has travelling wave solutions 



231, and 



The contraction of a cylinder end has been observed for the bacteria 
the travelling waves have been discussed in other contexts 0. Here we discuss the 
instability of the recoiling end and the final spacing of the spheres. 
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FIG. 14. Numerical simulation of a collapsing cylinder with a free end, showing contour 
plots of the density. As the cylinder recoils, aggregates are left behind. The spacing between 
the aggregates is determined by the density of the cylinder as it collapse. The final frame in 
this figure shows the density profile along the centerline of the cylinder for the last contour 
plot (with 4 aggregates). Note the maximum density is near the cutoff density = 500 
discussed in the text. 



Figure shows a set of numerical simulations of a retracting cylinder. The 
simulation shows that the "end" of the cylinder collapses into a spherical aggregate, 
and simultaneously, in front of the aggregate "waves" travel into the bulk of the 
cylinder. In the simulations, the ends of the cylinder actually try to collapse to 
an aggregate of infinite density. In order to continue beyond this singularity and 
simulate the formation of an array of aggregates we introduced a cutoff which emulates 
the biochemistry in the actual experiment: when the bacterial density becomes too 
high, the bacteria consume all the food sources (succinate and oxygen) in their local 
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environment and cease producing the attractant aspartate. We modelled this by 
changing the equation for c to 

at 

Here we have introduced p*, the cutoff bacterial density. In [0 this cutoff density 
was estimated (assuming the cutoff was caused by oxygen depletion) and shown to 
vary exponentially with the overhead oxygen concentration in the cell: p* ~ 6*"°^. 
For the simulation shown in Figure inithis cutoff is p* = 500. Note that the density 



of the (undisturbed) cylinder in front of the retracting rim slowly approaches the 
cutoff density p*; we have found in simulations that the undisturbed cylinder always 
collapses to a density close to the cutoff value. 

The formation of the density wave occurs because the retracting end perturbs the 
cylinder in front of it. We can find the time-evolution of perturbations to the cylinder, 
requiring that they decay away from the free end. The most unstable mode can be 
found using methods of stationary phase(cf. P8| ). 



If the linear growth rate is u;(g) the point of stationary phase satisfies 

duj , 
Im — =0 
dq 

duj Re(ti;) 

Re 



dq * Im(g) 



For the discussion here, we perform the calculation using the free space dispersion 
relation. A perturbation to constant density has the form 

p = Po + 5e-*-'?^ 
c = Pot + xe"*-''^ 

where we have taken e = 1 for simplicity in this calculation. Plugging into the 
equations and linearizing gives the dispersion relation: 

^ = -iy/poq + q^ 

The most unstable mode is, in dimensionless units, 

q* 
to* 
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These formulae demonstrate that the wavelength of the modulations are governed 
by the undisturbed density in front of the rim. Since this density is asymptotically 
determined by the cutoff p*, it follows that the wavelength of the ripples is deter- 
mined by the cutoff. The characteristic distance between the aggregates is, therefore, 
determined by the cutoff. This conclusion can be experimentally tested. 

The predictions for the wavelength and velocity of the front compare well with 
numerical simulations. On decreasing from 500 to 140 the wavelength of the ripples 
increases from ^ 1.4 to 2.5, in qualitative agreement with the formulas. 

We remark that the basic scenario outlined in this section was discovered by Elena 
Budrene, in unpublished experiments: After observing the travelling band to collapse 
as a cylinder (Fig. 1) Budrene observes fast "waves" propagating around the cylinder. 
Then, a fragmentation front moves across the collapsing cylinder, leaving spherical 
aggregates behind. The present theory predicts a scenario that is at least qualitatively 
very similar: the fast "waves" correspond to the excitations of the cylinder from the 
phase equation (i.e. the travelling steps, described in section 3). The fragmentation 
front occurs due to the mechanism outlined in this section. Unfortunately, it is not 
currently possible to make a quantitative comparison of the experiments to the present 
theory, though such a comparison would prove most interesting. 



V. CONNECTION TO EXPERIMENTS 



This paper has shown how the patterns formed by E. coli are connected to the 
geometry of singularity formation in the hydrodynamic description of the bacteria. 
We have constructed the solution for critical (two-dimensional) collapse, and devel- 
oped a theory for modulations to the cylinder. The phase equation provides a useful 
simplified description of a perturbed cylinder. Ultimately, the spacing of spherical 
aggregates is determined by the instability of a pinched cylinder of bacteria. 

Here we compare our work to published experiments and suggest tests of the 
theory. Not all the coefficients in the original equations have been precisely mea- 
sured 1^ for the experimental regime of interest. In particular, neither the attractant 
production rate a nor the chemotactic coefficient k have been measured for bacte- 
ria in the same chemical environment as that of the collapse experiments. Thus, at 
this stage we can make only order of magnitude numerical comparison with experi- 
ments. Here we use the values of the coefficients for bacteria in a liquid medium [0: 
bacterial diffusion coefficient = 7x 10~^cm^/sec; attractant diffusion coefficient 



Dj, = 10 cm /sec [49|; chemotactic coefficient k = 10 cm /sec; and attractant 



production rate a = 10 /second/bacteria 

An important prediction of this theory is the critical mass of bacteria for cylindri- 
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cal collapse: assuming the above parameters, the formation of a collapsing cylinder re- 
quires a minimum number of bacteria per unit length M = ADbDc/ (ka) = 3xl0^/cm. 
The existence of a critical number of bacteria for cylindrical collapse has been inferred 
from experiments [0] , but this number has never been directly measured. We empha- 
size that (with precise experimental measurements for the parameters Dh, k and a) 
the theory rigorously and precisely predicts this critical mass, allowing a direct test 
of the theory. 

In this paper we have extensively discussed how to construct the correct descrip- 
tion of two dimensional collapse. In the experiments, the subtle corrections to the 
dimensional scaling are probably not directly observable. However, the basic scal- 
ing relations expected from the similarity solution — for example, that the maximum 
density is related to the length scale of density variations by pm ~ L'"^ — could be 
measured in experiments, both for cylindrical and spherical collapse. So far, no quan- 
titative and controlled measurements of the bacterial density have been performed. 

In the section on cylindrical collapse, we discussed both the "transient" initial 
regime, and the final asymptotic regime which appears after true scale separation 
between the inner and outer solution is achieved. In practice, we expect that the ex- 
periments probe only the transient regime. The hard upper limit on bacterial density 
is when the bacteria are closely packed; this corresponds to a density of 2xl0^/cm^. 
A typical initial density of bacteria in the liquid experiments [0 is 10^/cm^. The 
density increases by at most three orders of magnitude during the collapse, in con- 
trast with the ten orders of magnitude it takes to reach the asymptotic regime in our 
simulations. 

To our knowledge, modulations to a collapsing cylinder have never been quanti- 
tatively measured in experiments (although as mentioned above, Budrene has made 
qualitative observations of this effect). The travelling and point singularities that 
we predict for a modulated cylinder may be observable. In particular, we predict 
that the radial and axial length scales should be different for the point singularity. 
Because the difference in these length scales arises directly from the corrections in 
two-dimensional collapse, a measurement of these length scales would test the validity 
of our derivation of the two-dimensional solution. 

The modulated cylinder ultimately pinches off a point. We have argued that the 
spacing of spherical aggregates is determined by the instability of a cylinder with an 
end. In practice, when does the modulated cylinder pinch off (forming an end)? To 
answer this question, we must know when our theory breaks down. Collapse to infinite 
density cannot happen for bacteria, because they have finite size. It was argued in 
that even before the hard packing density of bacteria is reached, oxygen depletion will 
stop bacterial collapse. Regardless of the specific mechanism for stopping the collapse, 
at some time the highest density part of the cylinder — the point singularity — will 
stop collapsing. This is the time of pinch off: the point singularity evolves much more 
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slowly than the neighboring, less dense regions of the cylinder. 

This argument about the point of pinch off gives a testable prediction of our model, 
because the spacing of aggregates depends on the maximum density of the cylinder. 
In dimensional units, the most unstable wavelength (and aggregate spacing) is 



where we have used values of the coefficients from above. Thus, varying the maxi- 
mum attainable density of the bacteria should cause the aggregate spacing to change 
according to this scaling law. In a formula for how the maximum density in a col- 
lapsed aggregate depends on the oxygen concentration Cqx was derived, and shown 
to be pm ~ 6*"°^. This implies that the wavelength of the pattern should decrease 
exponentially with the oxygen concentration; systematic experiments could test this 
prediction. 

The one solid prediction that can be compared with present experiments is that 
there is a lower bound on the aggregate spacing, that follows from the hard-packing 
density of bacteria. Using the characteristic size lO/xm of E. coli this is approximately 
10^/cm'^. Thus the measured aggregate spacing should always be above the lower 
bound 



where in this estimate we used the assumed values for the constants Dh, k and a. This 
lower bound agrees with experiments, in that the spacings are typically measured in 
millimeters. 
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APPENDIX A: REMARKS ON NUMERICAL METHODS 

The partial differential equations described in this paper were solved using second 
order in space, finite difference methods, supplemented with adaptive mesh refine- 
ment. The time discretization used a 6 weighted Crank-Nicolson-type scheme (i.e. in 
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the equation f — hf the right hand side is evaluated at time (n + 9) At, where At is 
the timestep). Typically, in the simulations with one spatial dimension, 6 = 0.6. For 
the simulations in two spatial dimensions, we used an ADI operator splitting method, 
which requires using 9 = 1. Because these methods are implicit, at each timestep 
a matrix inversion was necessary. This is the most expensive part of the numerical 
method. 

The most subtle aspect of the numerical simulations reported in this paper is the 
mesh refinement. Without good mesh refinement it is impossible to get close enough 
to the singularity to resolve the (logarithmic) corrects uncovered for the cylindrical 
collapse; without good mesh refinement in the two-dimensional simulations it would 
be impossible to acquire enough decades of data to test the phase equation-theory 
presented in section 4. The mesh refinement was performed as follows: 

One Spatial Dimension: The philosophy of mesh refinement employed in this paper 
(first explained to us by Jens Eggers) is to frequently implement gradual changes in 
the mesh, as opposed to infrequently implementing large changes. Mesh refinement 
is implemented every time the maximum density increases by one percent. During 
the refinement, the characteristic scale over which the solution varies is determined, 
and a mesh is constructed to well refine this scale; typically, this involves making 
sure there are about one hundred mesh points across the region that the solution 
varies significantly. The solution at the new mesh is constructed by cubic spline 
interpolation of the old mesh. Because the mesh is refined frequently, the changes 
to the solution occurring during refinement are small, and there are no convergence 
difficulties after refinement. The algorithm allows finding solutions over essentially 
arbitrary changes in the bacterial density with as little as few hundred mesh points. 
The simulations reported in the paper typically use a few thousand mesh points, in 
order to accurately compute the slow approach to the log log scaling regime. 

Two Spatial Dimensions: In two dimensions, the equations are solved using stan- 
dard operator splitting techniques. The mesh is rectangular, described by two func- 
tions, the X coordinate Xj and the y coordinate yj. Because the one dimensional 
algorithm described above can achieve arbitrary resolution with a few hundred mesh 
points, it is possible to well resolve density changes in both spatial directions using of 
order 10^ mesh points. The algorithm for these simulations works analogously to that 
for the one spatial dimension case described above: Every fifty time steps, the x grid 
(or y grid) is rcmeshed, in accordance with the criterea outlined above. Typically we 
stagger the remeshing between the two directions by twenty five timesteps. 

Both the one and two dimensional codes were tested extensively by checking the 
solutions against known analytical results. All of the results that are presented in 
this paper follow the philosophy that numerical results are only believable if they can 
be replicated by asymptotic solutions of the equations; in turn, asymptotic results 
are only useful if they show up in numerics. For the two-dimensional code, one might 
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worry that the operator sphtting coupled with the remeshing induces artifical biases 
in the numerics; besides checking our numerical solutions against solutions to the 
phase equation, we have also tested the two dimensional code by checking that it can 
reproduce the scalings and the similarity solution for spherically symmetric collapse, 
where the solution is known very well. 



APPENDIX B: PHASE EQUATION FOR A BACTERIAL CYLINDER 



In this section we fill in the details of the calculation of the phase equation. 
Evaluating the coefficients in equation (|20| ) we find ||5^ 



F{Ro, Co) = Ro + 
G{Ro, Co) = Ro + 



vK 

2 

vR'o vRqC'o 

2 2 
1 



H{Ro, Co) = 2Ro + ^i7vR'o + v'K ' ^vRoC'o - v'iRoC'oY) 



To evaluate the solvability condition, we need to find the zero mode of the adjoint 
to the linearized operator. In this case, the linear operator is the matrix 

V2 _ V ■ (VCo-) -V ■ (Ro-) 
1 V2 

All of the terms in A are self-adjoint except those of the form V ■ (VCq-). Under the 
definition (for a cylinder of bacteria) of the inner product (/, g) = J rdr J dzf*g we 
can determine the adjoint 

[V ■ (VCo-)]^ = -drCodr 

which gives the adjoint linear operator 




V^ + VrCoVr 1 \ 

-V-(/?oV-) v^J- 



This linear operator possesses a simple zero mode: A^(1,0) = 0. The coefficients of 
phase equation thus become 

ci = {l,F{Ro,Co)) 

C2={l,CiRo,Co)) 

cs = {l,H{Ro,Co)) 
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A subtlety comes when we evaluate the inner products: we must integrate (in 
similarity variables) to the upper limit of validity of the similarity solution. For a 
cylinder, this upper limit is ry*. the radius at which the solution matches onto the outer 
solution. From the asymptotics discussed earlier, we use that ~ A~^/^ — f{s). 
Evaluating the inner products, and taking the limit r — > 0, we arrive at the result 

-4 

C2 = 4 
C3 = -4 



Which gives as the phase equation 

7 



2 

I 
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